library(foreign)

dorm <- read.dta("dorm results.dta")

d <- NA

assign.dorm <- round(runif(1310, 1, 32))

sapply(1:32, function(x) sum(as.numeric(table(dorm$ethnic[dorm$house_num == x])/sum(table(dorm$ethnic[dorm$house_num == x])))^2))

d <- NA 
sd.d <- NA

for (i in 1:5000) {
  
  assign.dorm <- round(runif(1310, 1, 32))
  x <- sapply(1:32, function(x) sum(as.numeric(table(dorm$ethnic[assign.dorm == x])/sum(table(dorm$ethnic[assign.dorm == x])))^2))
  div <- x[dorm$house_num]
  div_01 <- (div - min(div, na.rm = T))/(max(div, na.rm=T)-min(div, na.rm = T))
  div_rev_01 <- 1 - div_01
  
  d[i] <- mean(div_rev_01)
  sd.d[i] <- sd(div_rev_01)
  
}

hist(d, xlab = "Diversity", main = "Diversity Score for \n 5000 Hypothetical Experiments", breaks = 50)
abline(v = .373, col = "red", lwd = 2)
